// Dolphin.cpp : Defines the entry point for the console application.
//

#include <iostream>

#include "ObjectiveFunctions/FEnergy.h"
#include "Includes/GlobalSearchers.h"
#include "Includes/LocalSearchers.h"
#include "Includes/Operators.h"

#include "Elements/Chromosome_Binary.h"
#include "Elements/Chromosome_Real.h"



#include "MemeticAlgorithms/MemeticAlgorithm.h"
#include "MemeticAlgorithms/MAMetaLamarckian.h"

#include "Rng/ihs.h"

using namespace std;

/// global variables
const char HELP[]="\n \
        NAME:  benchOF - benchmark forcefields \n\
        SYNOPSIS \n\
                benchOF -i[format] [input file]  [OPTION] \n\n\
        DIPSCRIPTION\n\
                [format]=  f ---> XYZ format \n\
                                        g03 ---> Gaussian output format \n\
                                        l ---> XYZ format of lattice \n\
                                        car ---> Accelys format file \n\
                                        \n\
                -ff [force field]: can be one of water forcefields \n\
                                OSS2 : (default) \n\
                                TTM2F, TTMT21F, TTM3F :  TTM family \n\
                                        or atomic models (openbabel required ) \n\
                                Ghemical,    //Ghemical force field. \n\
                                MMFF94, //   MMFF94 force field. \n\
                                MMFF94s,   // MMFF94s force field. \n\
                                UFF,    // Universal Force Field \n\
                                        or lattice forcefields (ttm_ice required) \n\
                                TTM21F_ICE, TTM3F_ICE, QSPCFW_ICE, SPCF_ICE \n\n\
                -a [parameter]:  parameter files \n\n\
                -u [unit]: unit=0,1,2 for HARTREE(default),KCALMOL and KJMOL, respectively.\n\n\
                -nBench: number of evaluations for benchmark (default = 1)\n\n\
                -nScale [nsize]: size for normalize, nBench will be nBench*(nsize/nAtom)^2 \n\
                -grad : benchmark gradient benchmark also (default is NOT) \n\n\
        NOTE \n\
                Version 1.0: using new OF and enable TTM2.1F and TTM3F \n\
                 \n\n";


void testMSLS(PotFitter* fit)
{
	int i, j;
	int nDim = fit->getNParam();
	int nDim_inp = fit->getNParam_inp();

	bool bound = false;
	int lslength = 20000;

	double* var = new double[nDim];
	double* fullvar = new double[nDim_inp];

	double* uBound = fit->getUBound();
	double* lBound = fit->getLBound();

	char stfile[100];
	char edfile[100];

	cout << "Start file: "; cin >> stfile;
	cout << "End file: "; cin >> edfile;

	ofstream f1(stfile);
	ofstream f2(edfile);

	cout << "Example: " << endl;


	fit->ConvertParam(fit->getParam(), fullvar);
	for(j=0; j<nDim_inp; j++) f1 << fullvar[j] << ",";
	f1 << fit->evaluate() << endl;

	cout << fit->evaluate() << " --> ";
	fit->RunLevMar(lslength, bound, 0);

	fit->ConvertParam(fit->getParam(), fullvar);
        for(j=0; j<nDim_inp; j++) f2 << fullvar[j] << ",";
        f2 << fit->evaluate() << endl;

	cout << fit->evaluate() << endl;

	int nSamples = 50;
	int* seed; *seed = time(NULL);

	int* points = new int[nSamples*nDim];

	cout << "Generating hypercube samples" << endl;
	ihs(nDim, nSamples, 5, seed, points);
	

	cout << "Multistart local search" << endl;
	for(i=0; i<nSamples; i++)
        {
                for(j=0; j<nDim; j++)
                {
                        var[j] = (uBound[j] - lBound[j]) / nSamples * ( -0.5 + points[i*nDim+j] ) + lBound[j];
                }
		fit->setParam(var);

		fit->ConvertParam(var, fullvar);
	        for(j=0; j<nDim_inp; j++) f1 << fullvar[j] << ",";
		f1 << fit->evaluate() << endl;

		cout << fit->evaluate() << " --> ";

		fit->RunLevMar(lslength, bound, 0);

		fit->ConvertParam(fit->getParam(), fullvar);
	        for(j=0; j<nDim_inp; j++) f2 << fullvar[j] << ",";
		f2 << fit->evaluate() << endl;
	
		cout << fit->evaluate() << endl;
	}

	f1.close();
	f2.close();
	cout << endl;
}

void testGA(PotFitter* fit)
{
	unsigned int i, j;

	ObjectiveFunction* f = new FEnergy(fit);


	f->nEvaluations = 0;

	Population<double> pop;

	double* params = fit->getParam();
	vector<double> indv;

	for(j=0; j<f->nDimensions(); j++) 
	{
		cout << params[j] << " ";
                indv.push_back(params[j]);
	}
	cout << endl;

	cout << "Test evaluation: " << (*f)(indv) << endl;

	Chromosome_Real* ch = new Chromosome_Real(f->nDimensions(), f->lowerBounds, f->upperBounds);
        ch->fromDoubleVector(indv);
	pop.push_back(ch);
	

	for(i=1; i<50; i++)
	{
		vector<double> indv;
		for(j=0; j<f->nDimensions(); j++)
		{
			indv.push_back(Rng::uni(f->lowerBounds[j], f->upperBounds[j]));
		}

		Chromosome_Real* ch = new Chromosome_Real(f->nDimensions(), f->lowerBounds, f->upperBounds);
		ch->fromDoubleVector(indv);
		
				
		pop.push_back(ch);		
	}
	
	Mutation<double>* mut = new Mutation_Gaussian(0.01);		
	Crossover<double>* crs = new Crossover_Uniform<double>(0.8);
	Selection<double>* sel = new Selection_RouletteWheel<double>(new Scaling_Linear());
	Recombination<double>* rec = new Recombination_KeepBest<double>(1);


	GeneticAlgorithm<double>* ga = new GeneticAlgorithm<double>(pop, f, mut, crs, sel, rec);
	
	//LocalSearch* ls = NULL;

	MemeticAlgorithm<double> ma(NULL);
	//MAMetaLamarckian<double> ma(ga);

	LocalSearch* ls1 = new LocalSearch_DSCG(f);
	LocalSearch* ls2 = new LocalSearch_DFP(f);
	ls1->stepLength = vector<double>(f->nDimensions(), 0.7);
//	ls2->stepLength = vector<double>(f->nDimensions(), 0.7);


	//ma.lsPool.push_back(ls1);
	//ma.lsPool.push_back(ls2);
	ma.ls = ls1;

	ma.ls->evaluationLimit = 10000;
	cout << "limit: " <<  ma.ls->evaluationLimit << endl;


	ma.pLS = 1.0;
	ma.maSelectionStrategy = ma.maLSBest;
	
	while(!ma.done())
	{					
		ma.evolve();
		//	if (ma.nGenerations() % 50 == 0)
		{
			cout << f->nEvaluations << " evals: " << f->bestEvaluation() << endl;
		}		

//		fit->RunLevMar(10000, false, 0);
	}

	cout << "Best solution so far: " << endl;
	for(unsigned int j=0; j<f->bestSolution().size(); j++) cout << f->bestSolution()[j] << " ";
	cout << endl << "Fitness: " << f->bestEvaluation() << endl;
}

int main(int argc, char* argv[])
{
	int nRuns=0;
	time_t duration=time(0);
	
	vector<string> arList;
	
	if(!isatty(0)){ string s; while(cin>>s) arList.push_back(s);}
	
	for(int i=1;i<argc;i++){ arList.push_back(string(argv[i]));}
	
	stringstream ss;
	for(int i=0;i<arList.size();i++){ ss<<arList[i]<<" "; }
	
	while(ss){
		string ar;  ss>>ar;			
		//parameters for working
		if(ar.compare("-nRuns")==0){ ss>>nRuns;}

		if(ar.compare("-h")==0){ cout<<HELP;return 0; }
	}


        //----------------

	PotFitter* fit = new PotFitter(arList);
        
	if (nRuns == 0) fit->evaluate(fit->getParam(),5);
	//else testGA(fit); //using levmar package to reparameterize
	else testMSLS(fit);

	duration=time(0)-duration;
	cout<<" Finished successfully in "<<duration<<" secs! "<<endl;

	return 0;
}

